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Liquid-liquid wetting failure is investigated in a two-dimensional Couette system with 
two immiscible fluids of arbitrary viscosity. The problem is solved exactly using a sharp 
interface treatment of hydrodynamics (lubrication theory) as a function of the capillary 
number, viscous ratio and separation of scale, i.e. slip length versus macroscopic scale 
of the system. The existence of critical velocities, above which no stationary solutions 
are found, is analyzed in detail in terms of the relevant parameters of the system. Com- 
parisons with existing analysis for other geometries are also carried out. A numerical 
method of analysis is also presented, based on diffuse interface models obtained from 
multiphase extensions of the lattice Boltzmann equation (LBE). Sharp interface and dif- 
fuse interface models are quantitatively compared face to face indicating the correct limit 
of applicability of the diffuse interface models. 



1. Introduction 

Despite many years o f resea r ch, the physi cs of moving contact lines ( De Gennesl ( 1985 ); 
Oron. Davis fc Bankofj (|l997l k lBlakel ([2006])) is still not completely understood. This lack 
of understanding stems from different factors. Dynamical wetting operates on scales ex- 
tending from the macroscopic to the molecular ones. In between these scales, with the 
usually small Reynolds numbers achi eved in such motion, the strong v i scous forces are 
balanced with surface tension effects ( De Gennesl (|l985l) : IVoinovl (|l976h : ICoxl (|l986l) L A 
dimcnsionless measure of this balance is provided by the capillary number Ca = [iU/o, 
comparing the viscous term at the contact line, r/U, with the surface tension, a, where rj 
denotes the liquid viscosity and U the contact line velocity. Liquid motion at finite capil- 
lary numbers induces changes in the shape of the interface, and the resulting macroscopic 
dynamic angle, Qm{Ccl), is different from its static equilibrium counterpart. Within this 
context, the main issue in contact line research is therefore to relate the macroscopic 
angle to the inner physi cs, w i th par t icular emph a sis on the mechanism s releasing the 
small scale sing ularities |Co2 (|l986h : iDe Gennesl (|l986h : IVoinovl (|l976h L In fact, it is 
well known that the viscous stress diverges at the contact li ne if s ome physica l mech - 
anism is not introduced to remove the singularity (Dussan ( 19791 ); De Gennesl ( 19861 ): 
Fismcn fc Pomeaul (l2000hV To overc o me this problem, numerous proposa l s cam e out 



Huh fc Scrivenl (ll97lMVoinovl |l976l): [Hocking fc Rivers! (|l982l) : lDe Gennesl (ll985l) : ICox 



19861 ): iDussan. Rame fc Garofl 



(|199ll )). all of them conducting to the introduction of 



a small length scale parameter £ s used to release the viscous singularity. An example is 



2 



M. Sbragaglia, K. Sugiyama and L. Biferale 



provided by the slip length at t he boundar ie s, where this small sca le parameter is related 
to the presence of a finite slip (Cox ( 19861 ): Huh fc Scriven ( 1971)). Anyhow, also othe r 
mechanisms su ch as intermolecular fo rces ( De GennesI ( 1986 ): |Pismen" fc Pomeaul(j2000l ); 
Jacaminl (|2000h : iDing fc Spelt! (120071) ) in the immediate vicinity of the contact line can 



be considered. 

Once the small scale singularity is removed, the macroscopic angle emerges as a func- 
tion of the capillary number 6m{Ccl): although the explicit forms may clearly depend 
on the specific model used to release the singularity, the common feature is that the 
sma ll capillary nu m ber l i mit s cales linearly with the contact line velocity, 6m{Co) ~ 
Ca (IVoinovl (jl976h : ICoxl (|l986l )). Differences are expected to emerge close to the wet- 



ting transition: when the liquid advances, a critical speed exists above which a sta- 



on the solid (Eggers 


(l2004al); 


Jacaminl ( 


20041): Snoeiier. Andreotti. Delon & Fermigier 


(2007|): Blal 


:e & Ruschakl (Il979l): Querel (1991): Podsrorski. F 


iesselles & Limat (200ll): 


De GennesI ( 


19861): Sedev & Petrovl ( 1991): Simnkins & Kuck ( 


20031)). The understand- 



ing of this transition is crucial. In fact, the breaking of stationarity can be interpreted 
as a lost of universality: the hydrodynamical regime does not support anymore a time- 
independent solution, indicating some singular behavior in the matching between inner 

and outer regions^ 

In his review, Kistler (Il993l) supports the assumption that wetting failure occurs when 
the dynamic angles reach zero degrees [J whereas the contact line is observed to become 
V-sha ped in the vicinity of the instability (jBlake fc Ruschakl ( 19791 ): Ghannum fc Esmaill 

di99i). 

Evidence for the existence of such critical points of cntrainmcnt has also been pro- 
vided by some recent theoretical works. The problem has been t ackled using a f ull hy - 
drodynamic ca lculation incorporating viscous effects on all scales (lEggerd (l2004al 120051) : 
Hocking! d2001)). In part i cular, within the framework of the 'Landau- Levich' problem 
(jLandau fc Levichl (|1942l ): iDeriaguin ( ^943h ). i.e. a plate plunging into or being with- 
drawn from a liquid bath, it has been shown ( Eggerd ( 2004aL 20051 )) that stationary so- 



lutions cease to exist above a critical capillary number Ca cr . The value of the emerging 
critical capillary number is quantified exactly: although we find universal features in terms 
of the microscopic angle (0 m ) of the liquid at the wall (Ca cr ~ 0f n ), a non- universality 
is also present in the geometric al prefac t or depending on the angle of inclination of the 
plate with respect to the liquid (jEggerd (|2004aL l2005h ) . 



The whole picture is also enriched by recent experimental observations (ISnoeiier. Delon. Fermigier fc Andreotti 
( 20061 )) where it is shown that the formation of solitary capillary waves can drasti- 
cally change the value of the critical speeds of entr ainment (pre-critical we t ting tr an- 
sition). Along these lines, linear stabil i ty an alysis ( Golestanian fc Raphaell ( 2001al lb): 
ISnoeiier. Andreotti. Delon fc Fermigier ( 20071 )) for the relaxation to external perturba- 
tions has also revealed that dispersion relations behave differently away and close to the 
critical point. 

Wetting failure has also been investigated in a recent paper bv lJacqmin ( 2004 ). where a 
liquid-liquid iso- viscous system in a Couette geometry was treated using various methods. 
This system, consisting of two walls moving with opposite velocities, has been treated 
in the p arallel flow approxima tion, with Fourier series method and also with phase field 
models ( Jacqmin ( 20001 20041 )). In these liquid-liquid systems, wetting failure has been 
found to set in well before the dynamic angle reaches zero degrees. 



f Throughout the paper notations are consistent with a macroscopic angle that decreases as 
the capillary number increases (see also the geometry in figure [T]with 9m smaller than 8 m ). 
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In this paper, we further explore and elaborat e on these issues by following a double- 
sided strategy. First, we extend the analysis of Jacqmin ( 2004 ) to arbitrary viscosity 
ratios, x = Ms/Wi where fi g and fii are a gas and liquid viscosity respectively. We then 
quantify the breaking of a stationary regime at a critical capillary number, Ca cr (x, A, 9 m ), 
which depends on the viscosity ratio, x, the microscopic wettability, 6 m , and the ratio 
between inner and outer scales, A = i s /H, where £ s is the slip length associated to a finite 
contact line slip and H is the distance between the two plates in the Couette geometry 
(see section[2]). The theoretical approach used is based on the lubrication approximation, 
i.e. a sharp interface treatment of hydrodynamics dealing with weakly bended interfaces 
( Oron. Davis fc Bankofl (1997])). In order to investigate macroscale geometry effects we 
also com pare ou r results to t hose of a similar analysis applied to the Landau-Levich 
problem (|Eggersl (|2004al . l2005| )). 

In the second part of the paper we will follow a computational pathway. We will study 
the sa me Couette flow with multiphase e xtensions of the lattice Boltzmann equation 
(LBE) (|Wolf Gladrowi (|2000h ; ISuccil (|200ll) ). In these diffuse inter f ace m et hods effective 



slip is induced by the finite width of the interface ( Seppecher ( 1996f k Ding fc Spelt 



(|2007M which is explicitly considered in the separation of two bulk phases. On one side, 
the LBE will serve as a benchmark of theory in those cases when the analytical approach 
is questionable (strongly bended interfaces and/or finite contact angles). On the other 
side, wc will use the sharp interface calculation to study the effects of finite width of 
the interface on the system. In fact, comparisons between a sharp interface theory (with 
slip-released singularity) and LBE (with diffuse-released singularity) is a way to probe 
the universality of the dynamics far from the contact line, i.e. that it does not depend 
on the microscopic details used to release the singularity. 

Dynamical benchmarks of LBE for multi- phase flows are also important, because 
of the difficulties to con trol corrections to the hydrodynamical limit in pres e nce o f 

strong dens i ty gradients (ISbragaglia. Benzi. Biferale. Sugivama. Succi fc Toschil (I2007TI; 

Shanl l2006l ): lLee fc Fischer! (l2006h : I Wagner! (l2003h : ICristea fc Sofoneal (l2003l ): lYuan fc Schaefer 



( 20061 )). We will show that LBE supports asymptotically the correct hydrodynamical be- 
havior, but it may be necessary to look on macroscopic scales much larger than the 
interface width. The parameter controlling deviations from sharp interface hydrodynam- 
ics is therefore going to be proportional to the ratio between the inner length scale (i.e. 
interface width) and the outer scale of the system. 



2. Lubrication Approach in the Couette flow 

The geometry we study consists of two parallel walls moving with opposite velocities 
ill w . The two walls have opposite wettability so that when the fluid is motionless the 
interface is a straight line angling from wall to wall (see figure[T]). For each fixed horizontal 
[j] location (x), the interface profile will be denoted with h{x). Computations will be 
carried out under the assumption of a finite and fixed slip length (^ s ), for simplicity the 
same on both walls. Other parameters are the capillary number, Ca = U w fjLi/a, and the 
viscosity ratio, x — fJ-r/fJ-i- It is further assumed that the left fluid is the more viscous 
one so that fii > fi r and consequently x < 1. The computation at x = 1 W1U serve as 
a benchmark t est for our results when compared to the results given in the paper of 



Jacqminl (|2004f ) 



Using a lubrication approximation ( Oron. Davis fc Bankofl (1997); Jacqminl ( 2004 ); 



f The origin of the coordinates is chosen so that the horizontal spreading of the interface is 
centered at x = 
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Figure 1. The geometry used to analyze the stability problem. Two parallel walls are moved 
with opposite velocities ±U W . The two fluids under consideration are distinguished in left (1, 
more viscous) and right (r, less viscous) and the viscosity ratio will be denoted by % = Mr/^i- 
The angles that the fluids are forming with the walls are complementary (opposite wetting 
properties). The microscopic angle 9 m is defined as the angle that the left fluid is forming with 
the lower wall and the angle 8m is taken as the angle in the center of the channel, with the 
convention that for finite capillary number 6m is smaller that m . The interface is determined 
by a function h(x) of the horizontal coordinate. Streamwise and vertical velocity fields for the 
left (right) fluid will be denoted with ui >x and W; >s/ (u r ,x and U r , v ). 



Eggerd (120051 l2004al) : ISnoeiierl (120051) 1 we carry out a long wavelength expansion in the 
stationary equations of motion (i.e. Stokes equation and continuity equation). Therefore, 
the theory developed is expected to be valid in the limit of small tilting angles between 
the interface and the wall. In this way, for a fixed horizontal coordinate x, the complete 
set of equations to be analyzed is 



/'/- 



d 2 u Lx 
dy 2 



dpi 
dx 



/'r 



dy 2 



dp r 
dx 



(2.1) 



d X Ul X + dyUly = d x 



d„u r .v = (2.2) 
and pressur e pi (p r ^ 



where we label with I (r) the left (right) fluid with viscosity /ij (fi 
These equations are the natural generalization of those analyzed in the paper of lJacqmin 
(|2004l ). where the author considered the case \ii — fi r . The solutions for the horizontal 
velocities in the left and right fluid u r _ x ) are rapidly evaluated as: 



1 ( A 1 
ui x = — [At + Biy + -pi, x y' 

W V 2 



f 

Mr 



1 

B rV + -=Pr,xV' 



(2.3) 



^ and p rt x 

with the usual boundary condition of zero normal velocity at 



T^p. The vertical components u^ y and u r ^ y have 



where we have used pi jX = 
to be evaluated from (|2.2 
the wall. Obviously, Ai, Bi,p^ x , A r , B r , p r ^ x have to be determined upon the imposition of 
ad- hoc boundary/matching conditions that will determine the position of the interface 
h(x). The relevant matching conditions are the continuity of the parallel velocity and 
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viscous stress at the interface 

i i 3u, : 



du r a 

h d v 



dy 

plus the lower and upper wall boundary conditions written as 

du r 



u, I -U dUl ' x 



(2.4) 



(2.5) 



H 



dy 

Finally, we have the kinetic condition of no flux across the interface 

fl X Ul tX + hyUl t y = h x U rtX + h y U rt y = (2.6) 

with n the normal at the interface in h(x). In Appendix A we show how to map the 6 
boundary conditions in a closed system and solve the corresponding matrix problem as a 
function of the separation of scale, A = i s /H, the capillary number, Ca — U w fii/j, and 
the viscosity ratio, \ = \i r j ' 

Once the pressure drop across the interface is known, one may derive using the Laplace 
law the equation for the local curvature, k, of the interface 

dn 

o—=Pr,x-Vi,x (2.7) 
dx 



with a the surface tension at the interface. As already noticed by I Jacamin ( 2004 ). if 



we define the interface arclcngth coordinate as s, we can write consistently with the 
lubrication approximation that 

dn 

0--J- = Pr,x ~ Pl,x (2-8) 

ds 

where the derivative of the curvature is connected to the angular variation ^ — n (where 
9 = 7r — 9). Summarizing, the governing equation set is: 

-^(K,9,x,y) = (o- _1 (p rx -pi iX ),-K, cos 6, -sin^). (2.9) 
ds 

The determination of the interface becomes now a non-linear boundary value problem: 
we need to solve the ODE (|2.9[) for a given capillary number, Ca, viscosity ratio, \i an d 
separation of scale, A = £ S /H, with the boundary conditions for the microscopic angle at 
the wall equal to 9 m 0. We have therefore a four-parameter problem. To solve numerically 
the previous set of non-linear ODE, we adopt a second order Rungc-Kutta method with a 
non uniform grid with increasing resolution near the wall. Notice that in presence of two 
different viscosities (/i r ^ /i/) the interface is not going to be symmetric with respect to 
the center of the channel. For a given value of viscous ratio x an d capillary number Ca, 
we look for solutions of (|2.9p by fixing the angle in the center of the cell, 9m, and choosing 
the curvature in that location so as to match with the desired boundary condition, 9 m . 

In figure [T] we show the results for 9m as a function of the capillary number Ca. 
The separation of scale is kept fixed to A = £ S /H = 10~ 5 and the viscosity ratio to 
X = 1-0. Various boundary conditions (i.e. microscopic wettabilities 9 m ) are considered. 
In the limit of small Ca the macroscopic angle is equal to the microscopic wettability 9 rn 
but, soon after Ca is increased, the interface is stretched and the macroscopic angle is 
decreasing. By increasing the capillary number, we clearly see that there is a range where 
two solutions can exist. This is typical of fixed point structures in dynamical systems, 



f for the present study we assume that the microscopic wall wettability is not dependent on 
the velocity 



6 



M. Sbragaglia, K. Sugiyama and L. Biferale 




1.0 

0.8 
§ 0.6 
^ 0.4 
0.2 
0.0 



. \ e M =3.i° y 

''•..(unstable) \ 


e*=25.4° 
/ (stable) 


A=10~ 5 
Ca=3.2- 10~ J 







0.002 0.003 
Ca 



0.004 




x/H 



Figure 2. Left: Macroscopic angle 9m as a function of the capillary number Ca for various 
microscopic angles 6 m and for a viscous ratio \ = 1-0- The value of the microscopic angle 
corresponding to each curve is provided by the small capillary number limit (Ca — » 0). The 
value of the separation of scale is kept constant to A = 10~ 5 . Data from Jacqmin (2004) are 
also reported for the case with 6 m — 45° (circles). On a specific curve (6 m = 45°), the critical 
capillary number, Ca cr , above which no stationary solution exists is indicated. Right: two static 
profiles in the stable and unstable branch for a viscous ratio x = 1-0. For the same microscopic 
contact angle (# m = 45°) and separation of scale (A = 10 -5 ) we show two stationary interface 
profiles corresponding to the same capillary number, Ca — 0.0032. The profile is plotted in 
terms of dimensionless coordinates, y(x)/H and x/H. Results correspond to the branch shown 
in figure[T]where our results are compared with those of Jacqmin (2004). Notice that the unstable 
solutions correspond to an enlarged bending of the interface close to the wall region. 



suggesting that a bifurcation happens at a critical capillary number, Ca cr , separating two 
branches, a stable branch < 0) and an unstable one (^ffj > 0). Above the critical 

capillary number no stationary solution can be found: beyond this value the interface 
has to evolve dynamically and liquid cntrainment on the solid is expected to take place. 
To better stress the existence of a range of capillary numbers where two solutions are 
found we also show (see right panel of figure [TJ the stationary interfaces corresponding 
to two different outer angles, 9 m, with the same boundary physics, 9 m . From this figure 
we also see t hat o ur analysis perfectly matches data presented in figure 2 of the paper 
bv I Jacqminl ( 20041 ) and reported in the left panel using the symbols. 

New results arc presented for the case of a different viscous ratio in figure [3] where 
we show 9m as a function of the capillary number for a given separation of scale, A = 
£ S /H = 10~ 5 , various microscopic wettabilities 8 m , and fixed viscous ratio, \ = 0-1. 
The behavior of the unstable branch is very sensitive to the viscous ratio as well as the 
microscopic wettability: here we see that as soon as x ^ 1, the second (unstable) branch 
docs not reach anymore very small values of 9m ■ In the right panel of the figure we show 
the critical capillary number as a function of the microscopic wettability. Notice that 
there are remarkable variations with respect to A only for the larger values considered: 
already for A = 10~ 5 — 10 -7 the results are pretty stable and weakly A dependent. 
This kind of study will be important later on, when comparing with the LBE where, 
because of numerical limitation, one cannot ever reach separation of scales smaller then 
A = 10~ 2 — 10 ~ 3 . Not surprisingly, figure [3] shows that for large scale separation the 
interface become more stable, i.e. Ca cr becomes larger. 

By changing the viscous ratio one can now understand the way Ca cr is related to x- In 
particular, in figurc|4]we present results for 9m as a function of Ca for a given separation 
of scale A = 10~ 5 at changing the viscous ratio, from x = 1 down to x — 10~ 3 (a realistic 
value for a liquid-gas interface at ordinary temperatures). The qualitative picture is 
always the same: the macroscopic angle decreases and reaches a bifurcation from which 




Figure 3. Left: Macroscopic angle 9m as a function of the capillary number Ca for various 
microscopic angles 6 m and for a viscous ratio x = O- 1 - The value of the microscopic angle 
corresponding to each curve is provided by the small capillary number limit (Ca — » 0). The value 
of the separation of scale is kept constant to A = 1CP 5 and the viscosity ratio to \ — 1-0. Right: 
the critical capillary number as a function of the microscopic angle 9 m for various separation of 
scales A: A = 10 _1 (times), A = 10~ 3 (squares), A = 10 -5 (circles), A = 10~ 7 (triangles). 



we extract the critical capillary number. We also notice that in each curve the range of 
macroscopic angles associated with the two branches is decreasing with x and also that, 
at decreasing Xi we reach a limiting curve, i.e. results become almost independent of 
the viscous ratio. The fact that by decreasing the viscous ratio we increase the critical 
capillary number can be qualitatively understood from the equations of motion (|2.1[) . 
When x ~~ > one of the two viscosities is getting close to zero whereas the other stays 
finite. By imposing a zero viscosity in the equations of motion we tend to rule out the 
contribution of the pressure gradient in the less viscous fluid. Therefore, the local change 
in the curvature given by equation (|2.8|) is triggered only by the pressure gradient of a 
single fluid, and we need higher capillary number to stretch it and break stationarity. 
All this is also complemented with the integral of the squared curvature along the profile 
(see right panel of figure 2]) for the same microscopic angle and separation of scale. The 
increasing stretching of the interface as a function of the viscous ratio (moving from 
X = to x = 1) is reflected in this plot. 
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Figure 4. Critical capillary number as a function of the viscous ratio. Left: for a fixed micro- 
scopic angle (0 m = 45°) and separation of scale (A = 1CP 5 ) we show the macroscopic angle 
8m as a function of the capillary number for different viscosity ratios X- Note the saturation of 
the curve while we are approaching values comparable with those of water and air at ordinary 
temperatures (x ~ 0.001). In that limit the equations are solvable exactly and we recover the 
usual solution for liquid-gas film. Right: We plot the integral of the squared curvature along the 
interface K = f \d xx h(x)\ 2 dx as a function of the capillary number for the same microscopic 
wettability 8 m — 45° and A = 10 _J . Different viscous ratios are chosen to stress the different 
interface stretching. 



2.1. Dependence on the walls wettability 

The aim of this section is to investigate the role of the wetting boundary condition at the 
walls. So far, we have considered the opposite wettability boundary condition as shown 
in figure [TJ where the angle the left liquid is forming with the lower boundary is the 
same as the angle the right fluid is forming with the upper boundary (configuration a). 
A possible test to understand the role of the boundary wettability is to consider another 
configuration (configuration b, see figure [5|) where the top wettability is kept fixed to a 
given value. The microscopic wettability is now changed only in the lower boundary. We 
can compute the critical capillary number in terms of 9 m and the separation of scale A in 
the same way as described in the previous sections. The separation of scales considered 
are between A = 10~ 5 and A = 10~ 3 . The critical capillary number in configuration b 
overestimates the case with opposite wettability (configuration a) and discrepancies are 
found to be enlarged when A is becoming larger (see figure [6]) . The separation between 
the two contact line regions is indeed less pronounced when A is not very small, and 
differences can be expected to emerge in that limit. Overall, the discrepancies between 
the two critical capillary numbers are not that large, but this is a consequence of the small 
values of the separation of scale chosen. Repeating the analysis for values of A ~ 10 _1 
leads to a larger discrepancy, indicating that non universalities with respect to the outer 
geometry become more and more pronounced by increasing A. 



2.2. The limit of zero viscous ratio 

We now address the properties in the limit of zero viscosity ratio, \ = Mr/w ~~ * 0- To 
do that we have to consider only the left liquid because the other component has zero 
viscosity and its prcfluctuations (p r . x = in equation (|2.ip ). In this limit the lubrica- 
tion equations simplif y and reduce to t he well known treatment for liquid-gas interfac es 
( Eggersl (|2005l . l2004al ): ISnoeiierl |2005l ): ISnoeiier. Andreotti. Delon fc Fermigierl (|2007t l). 
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Figure 5. Role of the geometry. The opposite wettability configuration (left, configuration a) is 
modified into another configuration (right) where the top wettability is kept fixed (configuration 
b). In both cases we take the macroscopic angle 9m as the angle corresponding to the location 
h{x) = H/2. 
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Figure 6. Outer scale (non) universality. The critical capillary number is computed for the 
case of an opposite wettability boundary condition (see configuration a, figure [5} and different 
separations of scales (Symbols). A similar analysis is done for configuration b (see figure^ where 
the top wettability has been fixed to a constant value, while changing only the lower wettability. 
The values of the separation of scales considered are between A = 10~ 5 and A = fO -3 and the 
viscous ratio is set to x = 0.1. 



where viscous friction is balanced with capillary effects: 

d 3 h _ 3Ca 
dx* ~ ~h(h + £ s )' 



(2.10) 



This is very simila r to the inner descriptio n found in a cont act line problem involving 

the Landau-Levich ( Landau fc Levich ( 1942);lDeriaguir] (Il943f)) geometry of a plate with- 

drawn at a given velocity from a liquid bath (lEggersI (l20051 ): ISnoeiier. Delon. Fermigier fc Andreotti 
( 2006( ) ; Snoeiier. Andreotti. Delon fc Fermigier ( 20071 )). This inner s o lution i s usually 



found by expanding i n a power series in the capillary number (jEgger 1 (|2004bl ): iHockingl 
(1983!); Voinov ( 19761 )) with the final result that, away from the contact line, the cube of 



10 



a macroscopic angle i (x) 
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3r„\ r~ I dh) 3 



, dx / 



is changing as 



\x) w^-9Colog 



(2.11) 



This is basically achieved for sufficiently large arguments of x/£ s and the related range 
is usually called intermediate range. When we are approaching Ca cr , the macroscopic 
angle is decreasing and one can argue that the relevant scaling properties of the critical 
capillary number can be captured by setting to zero the lhs of equation (|2.1ip . This 
implies that Ca cr is rescaling like the cube of the microscopic wettability and inversely 
proportional to log A -1 . Therefore, for a fixed microscopic wettability 9 m , it should be 
expected that by changing the separation of scale there is well precise rescaling factor 
for the critical capillary numbers 



Ca c 



Ca, 



cr | Ai 



log Ai 
log A 2 



(2.12) 



It is interesting to check the range of applicability of the previous heuristic argument 
with respect to our data. In particular, in figure [7] we check the validity of the proposed 
scaling relation for two characteristic viscous ratios, x = 1-0 and x = 0.1. We consider 
the critical capillary number as a function of the microscopic wettability 9 m and various 
small separation of scales A. We compare both unsealed critical capillary numbers and 
scaled ones (using equation 12. 12p . Interestingly enough, we observe that the scaling be- 
haviour is correct even for pretty large 9 m whereas the theory de veloped is base d on a 
lubrication approximation where small titling angles are assumed ( Snoeiier ( 20051 )). It is 
interesting to observe that the scaling relation is found to hold also in the limit of equal 
viscosities. The term proportio nal to log A i s connected to the viscous s tress divergence 
of the fluid at the contact line (jCoxl (|l986l ); |Pe Gennesl (jl986l ); IVoinovl (|l976l )) and this 
happens for all the values of the v iscou s ratio . Even if the equivalent of equation (|2. 1 1[) 
for finite viscous ratios is different ( Cox (1986)), the separation of scale is always appear- 
ing as log A and that determines the main rescaling properties for the critical capillary 
number. 

In the limit of x ^ it is also worthwhile to compare our estimate for the criti- 
cal capillary numb er in the Couette ce ll with a similar analysis done in some recent 
papers by Eggers ( Eggersl (l2004a . 20051)) fo r the case of the Landau-Levich geometry 
( Landau fc Levichl (1942); iDeriagufnl ( 19431 )) of a solid withdrawn from a liquid bath. 
In that case, the problem has been tackled as a multiscale problem: the inner contact 
line region where viscosity is balanced with capillarity must be connected with an outer 
region where gravity is balanced with capillary forces. An instability occurs when the 
outer meniscus approaches the shape corresponding to a perfectly wetting fluid, i.e. an 
apparent contact angle approaching zero degrees. An investigation of the conditions un- 
der which the highly curved contact line region can b e matched to the oute r profile leads 
to the introduction of a critical capillary number as (jEgger s1 (l2004aL[2005T) ) 



Ca c 



:i 

m 

9 



log 



Ca, 



1/3/ 



l&l^[M{s max )Y\9 v 



(2.13) 



where 9i n is the angle of inclination of the plate with respect to the liquid bath, Ai is 
the Airi function and s max = —1.0188... is the point where the Airi function assumes its 
maximum. The setup is clearly different from our Couette cell but neverthlcss is useful to 
compare both results to understand the role of the geometry in determining the critical 
capillary number. We choose the same separation of scale, A = 10~ 7 , and a small viscous 
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Figure 7. Scaling for the critical capillary number as obtained from the full hydrodynamic cal- 
culation in the lubrication approximation. Left: the critical capillary number Ca cr as a function 
of the microscopic angle 9 m for \ — 1.0 is rescaled according to the formula predicted by theory. 
A separation of scale is kept fixed to A2 = 10~ 7 (triangles). Various Ai are then considered: 
Ai = 10~ 3 (squares), Ai = 10~ 5 (circles). In the inset we show unsealed data (same symbols) for 
the critical capillary number Ca cr as a function of 6 m '- notice that the variability is extensively 
reduced for rescaled variables (i.e. good collapse). Right: Same as left figure with a viscosity 
ratio x = 0.1. 



ratio, x = 0.01, in the Couette cell. This is done to be in the region where we observe the 
collapse of data, as discussed in figure 0] We use two characteristic angles of inclination 
9 i n in (|2.13|) (9 i n = 5.73° and 9i„ = 90°) and the results are plotted in figure [8] The order 
of magnitude of the different critical capillary numbers is comparable and data share the 
same scaling properties with respect to 9 m (Ca cr ~ #f n ). Anyhow, the presence of an 
overall prefactor dependent on the geometry is clearly visible from the figure. For the 
sake of completeness, we have also considered a recent prediction of iHockina (|200l[ ) for 
fluids in narrow channels. The geometry is that of a fluid confined between two parallel 
plates at distance 2d under the effect of gravity. The main difference with the previous 
cases is that the speed of withdrawal is not an additional parameter but is fixed in terms 
of the bulk forcing, i.e. gravity. After introducing a Bond number, Bo = plB , with pi 
the liquid density and g the gravity, one can find the connection between Bo and the 
capillary number to be Ca = Bo/3. In this notation /xj is the liquid viscosity and U is the 
average speed of motion due to gravity. A critical value in the Bond number (equation 
(17) in Hocking(2001)) is then translated into a critical capillary number as 

The plot of (|2.14p is also considered in figure[5]with the same A of the previous estimates. 
This prediction turns out to overestimate all previous curves. 
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Figure 8. A comparison between our analysis and the proposed estimate for the critical capillary 
number given in equation (|2.13[) for a plate being withdrawn from a liquid at a given speed. 
The angle of inclination of the plate is 0i n and, as it emerges from Eggers (2004), the critical 
capillary number is a function of both the geometry and inner physics. Typical outcomes for 
din = 5.73° and 90° are shown. In our case we choose the same A and the viscous ratio is set 
to a small value so as to have data already collapsed on the limiting curve for \ — > (figure [4]) . 
Data from Hocking (2001) are also reported for the case of a fluid in a narrow channel under 
the effect of gravity. 



3. Numerical Approach: diffuse interface Methods and Lattice 
Boltzmann equation (LBE) 

In this section we will use numerical simulations to further investigate the problem pre- 
viously discussed. An approach based on first principles to address the inner physics of 
contact line problems is molecular dynamics (MD). These simulations, which usually in- 
volve Lennard- Jones liquids and sever al thousand molecules, appeared to exhibit contin- 

uum b e havior at the macroscopic lev el ( Thompson fc Robbing! (I 19891 ) ; Koplik. Banavar fe Willemsenl 
(19891); iDenniston fc Robbind (|200ll ): iBarrat fc Bocquetl (|l999l )). All the results of these 
molecular dynamics simulations offer insight especially in the region close to the con- 
tact line, pointing to the breakdown of the no-slip bou ndary condition at very sma ll 
distances as a poss i ble ex planation of contact line motion ( Thompson fc Robbinsl ( 19891 ): 
IBarrat fc Bocquetl (|l999h ) . Alternative explanations of contact line motion exist and 
they do not rely o n the breakdown of the no slip boundary condition. As was noticed by 
Seppecher ( 1996h . a sharp interface model may be questioned. It is itself an approxima- 



tion which may not be valid in the vicinity of the contact line. Then one may use molecular 
simulations or a continuum model able to describe the interface as a diffuse layer, i.e. 
a layer of finite thickness. The curvature of the interface near the contact line leads to 
mass transport across the i nterface, thus r e leasing the viscous singu l arity. Similar analy- 
sis in binary fluid systems (|jacqminl (|2000l ): IChen. Jasnow fc Vinalsl (|2000l )) showed that 
diffusive transport of the fluid leads to effective slip at the contact line. Also mesoscopic 
diffuse interface mod e ls for two-p hase flows based on the lattice Boltzmann equation 
( Wolf Gladro^(|200d) : ISu7dl (l200ll)1 were applied to t he sim u lation of contact line motio n 
and related problems (|b riant . Wagner fc Yeomansl (l2004h : iBriant fc Yeomand (2004 ): 
Latva-Kokko fc Rothmanl (|2007l ): ljia. McLaughlin fc Kontomarid (|2006 ): IZhang fc Kwokl 
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. n particular wc further elaborate along these lines using mesoscopic diffuse in- 
terface models where m ulti-phase phys i cs is induced by using a pseudo-potential approach 
originally developed bv lShan fc Chenl (|l993l . Il994l ). hereafter SC. 



3.1. The Lattice Boltzmann Equation for non ideal fluids 



We start from the usua l lattic e Boltzmann e quation with a single-time relaxation (jBhatnagar. Gross fc Krook 
(ll954h : IWolf Gladrowl (l2000h : ISuccil (120011) 1: 



f l (x + c l At,t + At)-f l (x,t) = -^(f l (x,t)-ti eq) (p,pu)') + F (3.1) 

where fi(x,t) is the kinetic probability density function associated with a mesoscopic 
velocity c; (there is a discrete set of velocities), r is a mean collision time (with At a 
time lapse), f[ eq \p,pu) the equilibrium distribution, corresponding to the Maxwellian 
distribution in the continuum limit and Fi represents a general forcing term whose role 
will be discussed later in the framework of inter- molecular interactions. From the kinetic 



distributions we can define macroscopic density and momentum fields as (jWolf Gladrow 
(|2000h : ISuccil (|200lh ): 



(3.2) 



For technical details and numerical simulations we shall refer to the nine-speed, two- 
dimensional 2DQ9 model ( Wolf Gladrow ( 2000l )), often used due to its numerical ro- 
bustness. The equilibrium distribution in the lattice Boltzmann equations is obtaine d 
via a l ow M ach number expansion of the continuum Maxwellian ( Wolf Gladrow ( 2000l ); 
Succil (1200 lb ) 



(eg) 



leg) 



C\pUi 



(44 



2cf 



-pu.Uj 



(3.3) 



where cj. — 1/3 and i = 1,2 = x, y runs over spatial dimensions. The we ights w^ eq ^ are 



chosen such as to en force isotropy up to fourth order tensor in the lattice (jWolf Gladrow 
( 2000f ); [Succil (|200ll )). Fro m the equilibrium dist ribution and the symmetry properties of 
q, it immediately follows ( Wolf Gladrow (^POOf O the kinetic second order tensor of the 
equilibrium distribution: 



where, in the first term of the rhs, we recognize the ideal-gas pressure tensor: Pjj = 
8ij(c 2 s p). In order to study non- ideal effects we need to supplement the previous descrip- 
tion with an inter-particle forcing. This is done by choosing a suitable Fi in (|3.1[) . In the 
original SC model, the bulk inter-particle interaction is proportional to a free parameter 
(the ratio of potential to thermal energy), Qb, entering the equation for the momentum 
balance: 

(3.4) 



= -Q b cl E w(h \ 2 )*P{x, t)i>(x + ciAt, t)cj 



being w(|c;| 2 ) the static weights and ip(x,t) = ip(p(x,t) the (pseudo) potential func- 
tion which describes the fluid-flu i d interactions triggered by inhomogeneities of the den- 
sity pr ofile (IWolf Gladrowl (|2000h : ISbragaglia. Benzi. Biferale. Sugivama. Succi fc Toschil 
(2007)). We shall refer to the pseudo-potential used in the original SC work, namely 



ip(x, t) = (l- exp(-p(x, t))). 



(3.5) 
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Note that this reduces to the correct form ip — > p in the limit p <C 1, whereas at high 
density (p 3> 1), it shows a saturation. This latter is crucial to prevent density collapse 
of the high-density phases (note that the SC potential is purely attractive, so that a 
mechanism stabilizing the high-density phase is mandatory to prevent density collapse). 
In principle, other functional forms may be investigated, sometimes with impress ive en- 
hancement of the density ratios supported by the model ( Yuan fc Schaefeil ( 20061 )). 
In order to understand the corrections to the ideal-state equation induced by the pseudo- 
potential, we need to define a consistent pressure tensor, P^, for the macroscopic vari- 
ables: 



-Fi + di(c 2 s p). 



(3.6) 



Upon Taylor expanding the forcing term and assuming hereafter At = 1, we obtain 



Ft 



(3.7) 



which is correctly translated ( Benzi. Biferale. Sbragaglia. Succi fc Toschi ( 20061 ): He fc Doolen 
(|2002l) ) into 



^Qbctdiipdjip. 



(3.8) 



The evolution scheme (|3.1[) together with the inter-particle force Q3.4p approximates the 
following diffuse interface equations for the density and momentum fields (|3 . 2[) : 



d t p + dj(puj) = 



(3.9) 



dt(pui) + Ujdj(pui) = -d 3 Pij + djUij (3.10) 

with I Lj j the usual viscous s t ress te nsor ( Wolf Gladrow ( 2000h : Succi ( 2001 ): Briant. Wagner fc Yeomansl 
( 2004 ): Briant fc Yeomansl ( 20041 )). Capillary effects, mimicking non-trivial interactions 

with the wall, can be incorporated in this appro a ch by using a suitable wall function 

Sbragaglia. Benzi. Biferale. Succi fc Toschi ( 20061 ): Sbragaglia. Benzi. Biferale. Sugiyama. Succi fc Toschi 



asa 
2007)) 



A disadvantage of the SC formulation is that, existing only one free parameter to tune 
(Ob), one cannot change independently the density ratio and the surface tension. Here we 

use a recent extension of the SC model that overcomes this limitation, by introducing first 

and se cond neighbors coupling in the pseudo-potential (jSbragaglia. Benzi. Biferale. Sugiyama. Succi fc Toschi 
(|2007t )): 



Fi = -c 2 s w(\ci\ 2 )Hx, t)[Giij>{x + c u t) + G 2 ip(x + 2ci,t)]4. 
i 

This coupling leads to the following Pressure tensor in the continuum limit: 



g 2 f iv^i 2 



^G2c%ipdj^. 



(3.11) 



(3.12) 



Furth er more, i n orde r to reduce the importa n ce of spur i ous cu r rents at interface s ( Sbra gaglia. Benzi. Biferale. Sugiyama. 

(120071): IShani (120061): iLee fc Fischerl (120061) : IWagnerl (|2003h : ICristea fc Sofoneal (120031) : 

Yuan fc SchaeferT ( 2006[) ) one can also introduce the non-i deal terms directly into the 

equilib r ium distribution as explain ed in some recent papers ( Briant. Wagner fc Yeomansl 

(|2004 ) : iBriant fc Yeomansl ( 2004 )). With this strategy we impose the desired pressure 

tensor instead of writing explicitly the forcing term (|3.11[) . We can thus use Q\ and Q 2 

to vary density ratio and width of the interface (i.e. surfac e tension) independently. As 

a matter of fact, the use of Van der Waals pressure tensors ( Briant. Wagner fc Yeomansl 
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( 2004 ; lBriant fe Yeomansl (2004|)), may limit seriously the density ratios achieved in the 



simulations. 

3.2. Stationary Contact Line description: Cox (1986) vs LBE 

In this subsection we explore and investigate the conditions under which our diffuse 
interface models converge to the predictions of sharp interface hydrodynamics. As we 
have already seen and discussed in equation (|2.11|) . we should expect an intermediate 
range with a prediction of scaling for the slope of the interface in terms of the capillary 
number and the microscopic slip length. A further extension fo r that predic tion of scaling 
for two fluids with finite viscous ratio has been provided by Cox ( 19861 ). This scaling 
relation translates (see figure 1 in Cox (1986)) in our coordinates (see figure [1} as 



9(0(r),x) = 9(0(0), X) + Ca g log(r/4) 



(3.13) 



9(0, X) 



f(x,x) 



dx 



(3.14) 



and 



f(x,x) 



2sinx [x 2 (x 2 



sin' 



x) + 2x~ 1 (x(tt - x) + sin 2 x) + ((-n - xf 



■ 2 
Sill X 



)] 



X 1 (x 2 



sin' 



x)((tt — x) + since cosx) + ((n — x)" 1 



x)(x 



where we have defined Ca„ 



sin x cos x) 
(3.15) 

as the capillary number estimated in the less viscous 



fluid, r as the distance from the stationary c onta c t line and the angle 0(r) consistently 
with the notation of figure 1 in the paper of lCol (|l986h . This is a pure sharp interface 
treatment of the moving contact line involving a finite viscous ratio and the length scale 
£ s is associated with microscopic slip motion. This turns out to be extremely useful to 
compare with our simulations where it is difficult to reach very small viscous ratios. 
We would expect that in our diffuse interface methods, when we look at the contact 
line on scales much larger than the interface width, £, we should recover a quantitative 
matching with the sharp interface prediction (|3.13[) . The length scale £ s is now associ- 
ated with the effective slip generated by the diffu se interface mechanisms ( Seppecher I 
(1996); I J acq mini ( l2000l) : iPismen fc Pomeaul ( 200(1) 1. This study opens the possibility to 
test the universality in the intermediate region independently of the inner physics mech- 
anisms. In our numerical simulations we fix the viscous ratio to be equal to x = 0.1 
and the associated density ratio is chosen so that the interface width £ is of the order of 
some grid points [j]. The wetting properties are chosen as explained in some recent studies 



Benzi. Biferale. Sbragaglia. Succi fc Toschi (120061) :! Sbragaglia. Benzi. Biferale. Sugivama. Succi fc Toschil 



2007t )). In particular, for the present set of simulations, we use a neutral wetting bound- 
ary condition (9(0) = ir/2 in (|3.13| ). To do that we assume that the pseudopotential 
of the lattice Boltzmann description (|3.5j) achieves a neutral stress at the boundaries 



(d n ip = 0, with n the normal with respect to the wall). With a fixed Ca g , interface width 
£ and viscous ratio x = 0.1 we then carry out numerical simulations in our Couette 
cell by enlarging the distance between the walls: H = 200Ax, 400Ax, 600Ax, 8OOA2; (see 
figure [9] left panel) , with Ax the lattice spacing. Half of the computational domain has 
then been initialized with the left, more viscous fluid and rotational boundary condi tions 
arc implemented at the inlet and outlet as discussed by iBriant fc Yeomansl ( 2004 ). No 
slip boundary conditions are implemented exactly, i.e. the mechanism of releasing the 



f One would need a systematic procedure to give £. A possible way is to define £ in terms of 
those points where the local gradient exceeds a given threshold. In any case the width £ does 
not exceeds some units and it is reassuring to assume (,/Ax — 4 ± 1 for the case of \ — 0.1. 
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Figure 9. The comparison between numerical simulations and the scaling relations for the 
intermediate region as predicted by Cox (1986). Left: The scaling for g(0,x) as a function of 
the distance from the contact line r. For a fixed capillary number in the gas phase, Ca g , fixed 
viscosity ratio, x> an d fixed slip length, the sharp interface theory predicts a linear scaling of the 
function g(0, x) (see equation (|3.13p ) with respect to log r . The development of the intermediate 
region is shown for various numerical resolutions approaching the sharp interface limit, where 
theory and numerics become asymptotically comparable. In the intermediate region the slope 
is the one predicted by theory, i.e. equal to Ca g — tJ,r ^ w . Right: the same as left figure with 
different capillary numbers and a given resolution in the vertical direction. 



contact line singularity is found in compressibility effects due to the diffuse nature of 
the model. The results of the numerical simulations (see figure [9]) clearly show that by 
enlarging the resolution, i.e. increasing the ratio of the outer scale H with respect to 
the interface width, we correctly approach the linear scaling behaviour of g(8(r), x) with 
respect to logr. The slope is correctly given by Ca g in (|3.13|) . The linearity with respect 
to Ca g is further checked in the right panel of figure [9] where wc keep fixed H and simply 
change the capillary number by changing the velocity at the wall U w . 

3.3. Full Hydrodynamic Calculations: Lubrication Theory vs LBE 
In the previous subsection we have quantitatively matched the diffuse interface meth- 
ods wrth_the universal scaling properties in the intermediate range, as predicted by 
Col (|l986l) . To achieve a good comparison, a lattice Boltzmann separation of scale, 



Albb = £/H, must be chosen as a small scale parameter, i.e. the scales of observation 
must be larger with respect to the inner length scale specified by the interface width £. 
Clearly, \lbe plays the role of the separation of scale in the diffuse interface descrip- 
tion. In particular, we can investigate the possibility of a quantitative matching with 
the prediction coming from the full hydrodynamical calculation developed in the first 
part of the paper. This would better explore the comparison between the two descrip- 
tions (sharp vs diffuse) on all length scales, and not simply in the intermediate range as 
done in the previous subsection. To do that we choose the same viscous ratio as before 
(x = 0.1) and we fix the distance between the walls to be H = lOOAx. The microscopic 
wettability is chosen to be equal to 9 m = 58° and then we vary the capillary number Ca. 
For each value of Ca we reach the stationary state and collect the corresponding value of 
the angle in the center of the cell (6m)- To extract the angle we do contour plots of the 
density field on the interface at the level p av = pg ^ p ' , with p 9t i the gas/liquid density 
respectively The corresponding plot of 6m vs Ca is shown in the left panel of figure [T0l 
Notice that when we reach a critical value of the capillary number, no stationary inter- 
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Figure 10. The comparison between theory and numerical simulations for the macroscopic angle 
9m as a function of the capillary number Ca on the stable branch. Left: for a fixed microscopic 
angle 9 = 58°, three different scale separations A are used: A = 0.048 (triangles), A = 0.024 
(circles), A = 0.012 (squares). In the numerics the different scale separations are obtained by 
fixing the interface width in the diffuse interface approach and using different heights H. Left: 
same as right figure for a microscopic contact angle 9 = 72°. 



face is observed above that value. This is precisely the same behaviour as predicted by 
theory. In particular, we can find the correct value of A in the sharp interface treatment 
of the Couette cell able to fit correctly the numerical prediction. It is found that the 
correct value of A is approximately A = 0.024, which is of the order of Xlbe- We can 
also repeat the whole experiment with the viscous ratio and interface width fixed to the 
same values, but doubling or reducing by a factor 2 the resolution in between the walls. 
We would expect and we observe (sec figure [10]) that the corresponding numerical data 
match with the theoretical prediction obtained by doubling or reducing by a factor 2 the 
previously used A. The same study is also done with a different microscopic wettability 
(@m = 72°, see right panel of figure [TO]) . Overall, we observe a good agreement between 
the two approaches. 

In figure fTTI we compare stationary interface profiles for a fixed viscous ratio, x = 0.1, 
and two characteristic capillary numbers, Ca = 0.022 and Ca = 0.033. We choose H = 
lOOAx and the interface width £ is kept fixed. The corresponding prediction from the 
theory has been produced with the value of A able to match the corresponding numerics 
in the plot of 9m as a function of Ca (see middle plots of figure fT0|) . It should be noticed 
that we match correctly the bulk regions, whereas small discrepancies emerge close to the 
boundaries. One may want to argue that the lubrication theory is valid only in the limit of 
small tilting angles (i.e. microscopic wettabilities), and the presence of finite wettabilities 
can produce a small mismatch between theory and numerics. Anyhow, it should also be 
noted that close to the boundaries the inner physics is completely different: in the sharp 
interface limit the dynamics is dominated by slip properties while in the diffuse interface 
case by the finite thickness of the interface. Overall, at least for the range of parameters 
studied, the agreement is satisfactory showing a good global univarsality with respect to 
the mechanisms resolving the contact line singularity. 

Much more stringent comparisons can be carried out looking at the local details in the 
stationary regime of the velocity field. Numerical and theoretical snapshots of the velocity 
vector are displayed in figure [T2j Qualitative features are very similar for both cases. In 
fact, it is evident from both theory and numerics the presence of a two- rolls structure 
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Figure 11. The static interface shape for a given capillary number and microscopic angle m - 
Left: for 9 m = 58° we show the results of the numerical simulations for two different capillary 
numbers: Ca = 0.033 (squares). Ca = 0.012 (circles). The corresponding theoretical prediction 
is shown (lines). The horizontal and vertical coordinate have been made dimensionless with 
respect to the system's height H . The separation of scale A in the theory is A = 0.024 that is the 
correct one to reproduce the diffuse interface results for the macroscopic angle 6m as a function 
of the capillary number Ca (see figure fTTJ|) . Right: same as left for m — 72° 



( Thompson fc Robbing! ( 1989| )) in the 2 fluids. Looking at the details of the velocity field 
in the streamwise (x) direction and we can further compare theory and numerics. In 
the left panel of figure [T51 for a fixed capillary number Ca, viscous ratio Xi separation 
of scale A and microscopic wettability 6 m , we show the streamwise component of the 
velocity (u x (x, y)) as a function of y for various horizontal locations x, as extracted from 
the theory developed in section [2] . Further details in terms of the separation of scales A 
arc displayed in the right panel of the same figure. When we fix the parameters Ca, x 
and 9 m , the velocity profiles tend to converge to the same master plot for small values 
of A. 

It is interesting to carry out a similar analysis in the numerics. To do that, we can fix 
the viscous ratio \, the microscopic wettability 8 m and the capillary number Ca to be the 
same as the case of the right panel of figurc[T2J With the interface width £ fixed to a given 
value, we can then increase the resolution between the walls and study the way the dif- 
fuse interface velocity converges to the sharp interface prediction. Results are presented 
in figure [Ml We notice that by increasing the resolution from H = lOOAx to H = 800Ax 
we reach already for H = 6OOA2; a limiting profile. On the other hand, if we compare the 
asymptotic numerical profile with the theoretical prediction (solid line in figure [T4|) we 
observe a mismatch in the less viscous region. The discrepancy can be due to different 
factors. First, we can argue that the sharp interface prediction coming from the lubri- 
cation theory may not be fully correct for finite wettabilities, i.e. a small overestimate 
in the pressure gradients would be amplified in the less viscous part as predicted by the 
stationary Stokes equations (|2.ip . It clearly would be possible to clar ify this point using 
a full 2-dimensional calculation, as for example the one discussed by ISomalinga fc Bose 
( 2000f ). Second, one should notice that in the interface region spurious effects emerge 
in the numerics. This is due to the presence of a stretched interface with loca l grad i- 
ents in the densi t y field (|Sbragaglia. Benzi. Biferale. Sugiyama. Succi fc Toschil (120071 ): 



Yuan fc Schaeferl (|2006l )). Neverthless, the comparison of the shape of the interfaces be- 



tween LBE and the sharp interface hydrodynamics did not reveal a great mismatch (see 
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Figure 12. A stationary state of the Couette cell in the numerics with diffuse interface (left) 
and in the sharp interface theory (right). The vertical length scale has been made dimensionless 
with respect to the height H and the horizontal one with respect to the horizontal extension 
of the interface L. Good agreement is found in the qualitative details of the velocity fields: a 
two-rolls structure is developed in both the gas and liquid phases. These plots share qualita- 
tive feature with the molecular dynamics simulations by Thompson and Robbins (1989). The 
interface location is also displayed (lines angling from wall to wall). 



x/L=0.25 
x/L=0.5 
x/L=0.75 





y/H y/H 

Figure 13. Velocity profiles in the streamwise direction as obtained from theory. Left: velocity 
field u x (x,y) as a function of y/H for different values of x along the horizontal direction (see also 
figure [T2"j) of the interface L: x/L = 0.25 (triangles), x/L = 0.5 (squares), x/L — 0.75 (circles). 
The separation of scale is A = 0.1, the capillary number is Ca = 0.011 and the microscopic angle 
is 9m = 58°. Right: The velocity profiles in the streamwise direction u x (x,y) as a function of 
y/H, for x/L = 0.5 and different separations of scales: A = 0.1 (triangles), A = 0.01 (squares), 
A = 0.001 (circles). 



figures 151 [TBI and ITTj) . i.e. to this regards the effect of spurious currents is certainly small. 
In order to progress in this issue, one should work out a better and refined numeri- 
cal scheme to reduce the spurious currents effect and check the c ompa ri son w i th the - 
ory again (Sbragaglia, Bcnzi, Bifcra le. Sugivama. Succi fc Toschil (120071): IShanl (l 2006 h : 
Lee fc Fischerl(|2006t) ; IWagnerl(|2003[ ): ICristea fc Sofoneal (|2003l ): lYuan fc Schaeferl pOOrj ). 
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Figure 14. Comparison between theory and numerical simulations. We show the velocity field 
in the streamwise direction u x (x,y) as a function of y/H, for x/L = 0.5 (L is the horizontal 
length scale of the interface, see also figure fT^)) . a capillary number Ca — 0.011 and a microscopic 
angle 9 m = 58°. The theoretical profile is obtained with a separation of scale A = 0.001 (already 
in the region of saturation shown in figure [T5|) . The numerical results are obtained with a fixed 
interface width £ and increasing numerical resolution in the vertical direction, i.e. decreasing 
the separation of scale \lbe = (,/ H. 



4. Discussions and Conclusions 

Wc have developed a sharp interface theory to describe a Couette cell consisting of two 
immiscible fluids. The two contact lines at the walls develop divergent viscous stresses 
and such a singularity is removed introducing a finite slip length at the boundaries (£ s )- 
The stationary properties of the system have been quantified in terms of the capillary 
number (Ca), the viscous ratio (x), the microscopic wettability (0 m ) and the separation 
of scale (A) between the inner physics (£ s ) and the outer geometry (i.e. the distance 
between the walls H). The problem can be closed using the assumption of small tilting 
angles at the interface (lubrication approximation), thus determining the stationary in- 
terfaces in terms of Ca, x, 9 m and A. It is observed that there it exists a critical capillary 
number, Ca cr , above which no stationary solution can be found. This critical capillary 
number corresponds to the case where liquid deposition occurs on the solid: the interface 
is stretched at the point that it cannot sustain any longer a stationary profile that is 
broken in favor of liquid cntrainment on the wall. Our analysis offers the possibility to 
examine this critical capillary numbers in terms of all the parameters, e specially t he vis - 
cous ratio, thus extending a similar analysis done in a recent paper bv ljacqminl (2004). 



Moreover, in the limit of s mall viscous rat ios, it allows to compare our results with the 
recent analysis proposed bv lEggerd (|2004af) for the case of a plate withdrawn at a given 



speed from a liquid bath. This elucidates the role of completely different geometries in 
determining the critical capillary numbers. 

In the second part of the paper we have studied the same system with a diffuse interface 
model based on the Lattice Boltzmann equation (LBE). The comparison with LBE re- 
sults allows us to benchmark the numerical model and to understand the effects of finite 
thickness of the interface on global quantities, as the critical capillary number, as well 
as local ones, as the interface shape and the velocity profiles. Good agreement is found 
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when the scale separation in the LBE, given by the ratio between interface thickness £ 
and distance between the walls H, becomes small enough. 

Extensions of the LBE simulations to other geometries, as it is the case for the Landau- 

Levic h case in pres e nce of gravity, would allow for direct comparison with experimental re- 

suits (lQuerel ( 1991 ): Snoeiier, Delon. Fermigier fc Andreottil (|2006h : ISnoeiier. Andreotti. Delon fc Fermigier 
( 20071 )) and also offer some more physical insight in problems that are still not com- 



1991 



pletely understood, as the importance of roughne ss, con t act angle hysteresis and speed 
depen d ency of the microscopic w e tting p roperties (iQuerel (Il99ll) ; Golestanian fc Raphael 
(20031): iRame. Garoff fc Willsonl (120041): iHeine. Grest fc Webbl &2004)). 

Very r ecently, a withdrawal experiment by Snoeijer et al. (jSnoeiier. Delon. Fermigier fc Andreotti 
(2006)) has been performed and, contrary to previous experiments ( Sedev fc Petrovl 



)), the transition to film entrainment did not occur at the critical capillary num- 
ber predicted by the theory (jEgger 1 (|2004aL [ioblh v This avoided critical bahaviour is 
due to the formation of a capillary ridge that does not trivially match the liquid film 
and determine the critical speed of entrainment. This can be related to contact angle 
hysteresis that has not been treated in the continuum models: the sensitivity of Ca cr 
with respect to 9 m (see also our figure [3]) can also support this interpretation. Another 
interesting possibility would be the introduction of a mesoscopic roughness in the nu- 
merical simulations and st udy the way this critical behaviou r changes in terms of the 
underlying heterogeneities (jKusumaatmaia fc Yeomansl (|2007j) ). One has to caution that 
such a kind of analysis would be limited to finite viscous ratios, unless one is opting for 
more sophisticated schemes incor porating density ratios comparable with li q uid-gas in- 
terfac es used in the experiments ( Inamuro. Ogata. Taiima fc Konishi ( 2004 ): Lee fc Lin 
(|2005h L Moreover, the properties of the diffuse interface models will emerge as a function 
of the separation of scale Xlbe- reaching extremely small values of Xlbe is somehow 
prohibitive because of the extremely large resolution needed to simulate a set of scales 
ranging from nm to tens of fim. To this regard, similar analysis to that presented in this 
paper (see section |2~2| would help to translate the numerical observations into realistic 
numbers. 



We are indebted with B. Andreotti, R. Benzi, J. Yeomans, H. Kusumaatmaja, F. Toschi 
and S. Succi for useful and enlightening discussions. 



Appendix A. Matrix problem for the lubrication approximation 

In this appendix we describe in details the calculations needed to solve the Couette 
problem in the lubrication approximation. Consistently with the assumption of small 
tilting angles, the normal vector to the interface is 

-^tJ (A1) 



with e x and e y Cartesian unit vectors. A zero normal component for the left field (for 
the right one similar arguments can be applied) at the interface means 

ui, n \ h = u y .i\ h - ^ u x ,i\ h = 0. (A 2) 

Moreover from the continuity equation (|2.2[) we derive that 

d [ h{x) dh 
Uy ^ h = ~dil u x> idy+ — u x ,i\ h . (A3) 
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The use of equation (|A 3j) together with equation (|A 2j) leads to 

d f h ^ 

Ul ^ h= dx J Ux - ld V- ( A4 ) 

Since in the inner (outer) region of the system there is no net mass flow rate, this means 
that 

r h(x) rH 

u x ,idy = / u x , r dy = 0. (A 5) 

Jh{x) 

The 6 boundary conditions that we need to fix Ai,Bi,pi^ x ,A r ,B r ,pi_ x are then represented 
by (|2.4p . (|2.5p and (|A5|) . They translate in the following system: 

h 2 ( h 2 \ 

A t + hBi + —p hx -x 1 [A r + hB r + —Vr,x J =0 (A 6) 

+ hp hx - B r - hp r>x = (A 7) 

Ai-£ s Bi=fnU w (A 8) 

A r + (H + t s )B r + + Pr,x = -fJ> r U w (A 9) 

h 2 h 3 

hA x + —Bi + —p ltX = (A 10) 

(H - h)A r + X -{B 2 - h 2 )B r + ^(H 3 - h 3 ) Pr , x = 0. (All) 
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